Extreme value Theory provides early warnings for critical transitions 
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Early warnings of critical transitions have been extensively used to detect abrupt changes of dy- 
namical regimes. In this paper we introduce new indicators based on the analysis of parametrically- 
modulated modifications in the properties of extremes for chaotic systems which possess experimen- 
tally accessible and detectable extreme value laws when far from bifurcation points. By measuring 
the deviation from the theoretically expected asymptotic distributions for long but finite samples, 
it is possible to detect the approaching critical transitions. Moreover, relations and connections 
between traditional and extreme value based indicators are explained and commented in detail. 
Numerical experiments have been performed on a stochastic differential equation describing the 
motion of a particle in an asymmetric double well potential under the effect of white noise. The 
results provide a gateway for using operatively the method described in other systems or data series 
analysis. 



I. INTRODUCTION 

An astonishing variety of complex systems ranging 
from finance to climate can experience abrupt changes 
at which a sudden shift of dynamical regime occurs. 
These so called "critical transitions" or "tipping points" 
can be formally seen as bifurcations in the dynamical 
systems terminology pQ. There is currently a great 
interest in understanding the dynamical behaviour in 
the proximity of a tipping point mainly for its impor- 
tance in analysing and forecasting events of social and 
economic relevance [2 A growing interest in this 
topic has emerged and many indicators of criticality 
have been developed in order to identify early warnings 
of abrupt transitions to different dynamical states. The 
most used indicators are based on the modifications of 
the auto-correlation properties of particular observables 
when the system is pushed towards a transition, usually 
accompanied by an increase of the variance and the 
skewness of the distribution of the observable analysed 
[2H1]- Although for specific systems such techniques 
successfully highlight approaching tipping points, several 
difficulties prevent to extend these results to general 
cases. In particular, for high dimensional and complex 
systems featuring oscillations and intricate bifurcation 
patterns, early warning signals may be misleading. The 
main issue for these systems concerns the accessibility 
of dynamical and geometrical properties of the physical 
measure (i.e. the attractor) in the proximity of tipping 
points and the difficulties in recognizing the nature 
of the bifurcation involved. Moreover, since all the 
early warnings indicators rely on the knowledge of 
asymptotical statistical properties of the systems, so 
called "false alarms" - wrongly identified bifurcations 
- may arise when considering finite dataset and model 
outputs [3]. 

In this context, it is highly desirable to build early 
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warnings indicators which provide themselves some 
estimations of dynamical and geometrical properties of 
the system so that false alarms can be discriminated 
from real early warnings. The main achievement of 
this paper is to introduce a new indicator of criticality 
based on the extreme value analysis which features the 
possibility to explore the asymptotic properties of the 
system and connect the fluctuations of maxima and 
minima of physical observables to the ongoing critical 
transition. To this purpose we exploit the theoretical 
results obtained so far for extreme values of dynami- 
cal systems to study the behaviour of a prototypical 
stochastic model in the proximity of a tipping point, 
showing analytically and numerically the capabilities of 
the new indicators. 

Although our aim is to provide a general tool to analyse 
tipping points, the motivation for the research presented 
here originates mainly from the analysis of the data 
series of turbulent energy in the so-called plane Couette 
flow, i.e. the flow of a viscous fluid between two parallel 
plates, one of which is moving relative to the other [5]. 
In this system the control parameter is the Reynolds 
number whose value determines either turbulent or 
laminar flow regimes. We observed that, when the 
Reynolds number approaches the threshold at which 
the turbulence decays, the probability of observing a 
very low minimum of turbulent energy increases. As 
a consequence, maxima and minima exhibit an asym- 
metric behaviour that can be quantified by applying 
the tool provided by the Extreme Value Theory (EVT). 
This experimental result has been the starting point to 
devise a heuristic approach for more general scenarios, 
which is the subject of this paper. The results on the 
ongoing analysis of the breakdown of turbulence in the 
plane Couette flow, performed in collaboration with P. 
Manneville, will be reported in forthcoming papers. 

The distribution of the largest or smallest values of cer- 
tain observables has been widely studied as it is of great 
interest in many practical applications (e.g. the analysis 
of environmental extreme events). This interest has led 
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to a fast development of a so-called extreme events the- 
ory [6j-[9]. The traditional approach introduced by Gne- 
denko [5] is based on the analysis of extremes obtained 
dividing the considered dataset into bins of fixed length 
and choosing the maximum for each bin. There, statis- 
tical inference is performed on such a reduced dataset 
by considering as model the Generalized Extreme Value 
(GEV) distribution family, which includes, as members, 
the Weibull, Gumbel and Frechet distributions, which 
greatly differ in terms of mathematical properties. Be- 
sides the block maxima approach, the so called Peak- 
over-Threshold (POT) is also widely used to tackle the 
problem of extremes. Extreme values are selected as ex- 
ceedances over a certain threshold, fitted to the Gener- 
alized Pareto (GP) distribution. The asymptotic con- 
vergence to the GP model is then guaranteed by the 
Pickands-Balkema-de Haan theorem [10] . 

Whereas the theory has been originally designed for 
the study of extremes for series of random variables, 
in the last decade the existence of asymptotic laws 
has been proven for maxima of observables computed 
on the orbits of dynamical systems. From the first 
rigourous paper by Collet in the period of a few 

years relevant new results have been obtained [T2"HT7] . 
The starting point of all these investigations has been 
associating to the stationary stochastic process given 
by the dynamical system a new stationary independent 
sequence distributed according to the GEV model, 
and then pulling back the obtained statistical laws to 
the original sequence extracted from the dynamical 
system. The assumptions necessary to observe a GEV 
distribution in dynamical systems rely on the choice 
of suitable observables and the fulfilment of particular 
mixing conditions. These results can also be used to 
study extremes of stochastically perturbed dynamical 
system as recently shown in [TH] [19] . 

In order to construct robust indicators of criticality, 
able to discern when the underlying system is close to 
a bifurcation, we will take two complementary points of 
view: 

1 . The classical approach of looking at the autocorre- 
lation properties of the time series of the considered 
variable is replaced with the investigation of the ap- 
propriateness of GEV for the extremes of the time 
series: the presence of long time-correlations is in 
antithesis with the fulfilment of the mixing condi- 
tions referred to above. 

2. The search for changes in the skewness of the 
bulk statistics is substituted with the analysis of 
whether parametric modulations to the system lead 
to changes in the qualitative properties of the ex- 
tremes, i.e. whether we observe transitions among 
the three kinds of distributions included in the 
GEV family. 

The capabilities of these new indicators will be anal- 
ysed for probably the simplest conceivable dynamical 



model featuring multi-stability, i.e. the overdamped ID 
motion of a stochastically forced particle in a double well 
potential. The choice of this model naturally relies on 
the fact that, even for complex and chaotic dynamical 
systems which feature tipping points, a reduction to 
the Langevin model is often attempted as it opens 
the way to approaching the problem in terms of the 
one-dimensional Fokker-Planck equation [30] , even if one 
must handle with care the construction of a surrogate, 
stochastic dynamics, as discussed in [21] regarding a 
relevant example of geophysical relevance. Moreover, 
this model can be easily explored numerically and 
analytically. In future papers we will test the indicators 
in complex systems arising from fluid dynamics. 

The fundamental reason why we believe that EVT- 
based indicators are preferable to usual methods based on 
the bulk of the distribution lies on the fact that, whereas 
there is, obviously, no universal model for the bulk statis- 
tical properties of a given system, such universality, even 
if with limitations, exists exactly for extremes for which 
EVT provides a universal family of parametric distribu- 
tions. 

This work is organised as follows: in Section 2 we 
present some basic notions of extreme value theory for 
independent and identically distributed variables and dy- 
namical systems. In Section 3 we introduce the new indi- 
cators and describe their properties. Section 4 is devoted 
to explaining the theoretical and numerical results for 
the motion of a particle in an asymmetric double well 
potential under the effect of noise. Finally in Section 5 
we review the results obtained for the system analysed 
and propose an algorithm to extend the analysis to other 
dynamical systems or time series drawing our conclusions 
and more suggestions for further work. 



II. ELEMENTS OF EXTREME VALUE 
THEORY FOR DYNAMICAL SYSTEMS 

A. Traditional Extreme Value Theory 

Gnedenko [5] studied the convergence of maxima of 
i.i.d. variables 

Xq, Xi, ..A m _i 

with cumulative distribution (cdf ) of the form: 

F m (x) = P{a m (M m - b m ) < x} (1) 

Where a m and b m are normalizing sequences and M m = 
max{V , Xi, X m ^i}. Under general hypothesis on 
the nature of the parent distribution of data, it has been 
shown that the distribution of maxima, up to an affine 
change of variables, converges in the limit for m — > oo 
to a single family of generalized distribution called GEV 
distribution with probability density function: 
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which holds for 1 + k(x — pi) /a > 0, using /jgl (loca- 
tion parameter) and a > (scale parameter) as scaling 
constants in place of b m , and a m [7] [IS], k € R is the 
shape parameter, also called the tail index, and discrim- 
inate among the classical Extreme Value Laws (EVLs): 
when k — y 0, the distribution corresponds to a Gumbcl 
type (type 1 distribution). When the index is positive, 
it corresponds to a Frechet (type 2 distribution); when 
the index is negative, it corresponds to a Weibull (type 3 
distribution). The type of distribution observed is very 
important as it discriminates the kind of tail decay of the 
parent distribution. A Gumbel distribution is typically 
observed when the parent distribution features an expo- 
nential tail, whereas Frechet and Weibull laws are typical 
associated to power law tails: a Weibull is usually ob- 
served when the tail is limited above by a certain thresh- 
old whereas a Frechet when the tail is limited below. In 
order to compare the properties of maxima and minima 
distributions, they should both be treated as they corre- 
spond to right or left tails of the parent distribution. In 
the forthcoming analysis, we will always change sign to 
the minima distribution |22) . 

As shown in Eq. [l] Gnedenko's results are related to 
a precise way of selecting extremes, the so called block 
maxima approach: it consists of dividing the data series 
of length s of some observable into n bins each contain- 
ing the same number m of observations, and selecting the 
maximum (or the minimum) value in each of them [2"2"] . 
The extremes are then fitted to the GEV distribution. 
In Section 4. A we will give a detailed description of the 
inference method used to extract the parameters of the 
EVLs. 



B. Extremes of Dynamical Observables 

We now consider the generalization of EVT to dynam- 
ical systems. The results presented in this section are 
provided in terms of maps (discrete-time dynamical sys- 
tems) but they are also valid for flows (continuous-time 
dynamical systems). 

Let us consider a dynamical system (f2, B, v, /), where f2 
is the invariant set in some manifold, usually M. d , B is the 
Borel cr-algebra, / : f2 — >• Q is a measurable map and v a 
probability /-invariant Borel measure. 
In order to adapt the extreme value theory to dynamical 
systems, we will consider the process Xq, X\, ... given by: 



where 'dist' is a distance on the ambient space Q, £ is 
a given point and g is an observable function. As we said 
above, we will use three kinds of observables gi,i = 1,2, 3, 
defined as: 



Si(a;) = -Iog(r) g 2 (x)=r- l> g 3 (x) = C - (5) 

where r — dist(x, £), C is a constant and (3 > E M. 
Using these observables we can obtain convergence to 
a distribution of type 1,2,3 if one can prove two suffi- 
cient conditions called D2 and D' that the dynamical 
system obeys. These conditions require the presence of a 
sufficiently fast decay of correlation for the stochastic dy- 
namical sequence and limit the possibility of having clus- 
tered extremes. Another way to approach the problem 
of extremes for the gi observables relies on studying the 
statistics of first return and hitting times, which provide 
information on how fast the point starting from a certain 
initial conditions returns in its neighbourhood, see the 
papers by |13j and |23| . They showed in particular that 
for dynamical systems preserving an absolutely continu- 
ous invariant measure or a singular continuous invariant 
measure v, the existence of an exponential hitting time 
statistics on balls around v almost any point £ implies 
the existence of extreme value laws for one of the observ- 
ables of type gi,i — 1,2,3 described above. The converse 
is also true, namely if we have an extreme value law which 
applies to the observables of type gi, i — 1, 2, 3 achieving 
a maximum at £, then we have exponential hitting time 
statistics to balls with center £. Recently these results 
have been generalized to local returns around balls cen- 
tred at periodic points [T5] . 

Since it is difficult to check the exponential decay of 
hitting time statistics as recurrences are hard to tackle 
analytically or numerically in more than three dimen- 
sions, it is highly desirable to connect the theory ex- 
plained above to a more straightforward dynamical in- 
dicator computable for a wide class of observables as 
well as in high dimensional dynamical systems. On the 
other hand, condition D2, introduced in its actual form 
by Freitas-Freitas [12] , could be checked directly by esti- 
mating the rate of decay of correlations for Holder observ- 
ables. Starting from this observation, Aytac et al. proved 
the existence of EVLs for dynamical systems whose cor- 
relations decay is at least summable [18] . Since these de- 
cays can be easily computed numerically, in the following 
discussion, instead of checking on the mixing conditions 
D2 and D' , we will refer exclusively to the results de- 
scribed in [T8] . 

In [22] and [25] the authors have verified the convergence 
to the classical EVL and the relations with the GEV fam- 
ily from both an analytical and numerical point of view 
in a wide class of mixing low dimensional maps, show- 
ing that - in mixing maps - GEV shape parameters only 
depend on the local dimension of the attractor 5((): 
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In [35] it is shown that a GEV distribution can be 
fitted only if the mixing conditions are fulfilled whereas 
other kind of distributions not belonging to the GEV 
family are observed for quasi-periodic and periodic mo- 
tions. The theory has been extended to stochastically 
perturbed dynamical systems in |18j : in particular, sys- 
tems whose trajectories under the deterministic dynam- 
ics evolve towards a finite number of fixed points possess 
EVLs when perturbed with additive noise. In these cases, 
the extreme value parameters depend on the phase space 
dimension and not on the local dimension as the noise 
helps the perturbed system to explore the neighbour of 
a fixed point in a ball which scales exactly with the di- 
mension of the phase space. In [TJj], the authors define 
some guidelines to effectively observe such EVLs in nu- 
merical simulations, pointing out the range of noise to be 
used. These results will be particularly useful for describ- 
ing the behaviour of the stochastic differential equation 
whose deterministic part is driven by a double well po- 
tential dynamics which will be extensively analysed in 
Section 4. 

In the next section we will describe what we expect to 
happen in a general case. The notations will refer to a de- 
terministic dynamical system which features at least two 
disjoint attractors but - under the conditions previously 
described - the results are valid also in the stochastically 
perturbed case. 

III. TIPPING POINTS AND EXTREMES 

Let us assume that v is the physical measure of 
/ : Q — > £1. In fact, since a dynamical systems may have 
several different invariant measures, we will always refer 
to the physical measure as the one arising naturally in 
any numerical simulations. We will also assume that 
for a given initial condition Xq, the physical measure 
v(xo) is unique but it is not the same for all the initial 
conditions in the phase space, i.e. there are at least 
two disjoint attractors. Let us perturb the dynamics by 
introducing a parameter A such that / = f(x ,X) and 
v = v(xq, A). We define A c as the critical value of A such 
that for < A < A c all the disjoint attractors continue 
to exist - each one with its own basin of attraction - but, 
whenever A > A c , the system undergoes a bifurcation 
that makes one of the disjoint attractors disappear. For 
/(A < A c ), the system has a summable decay of correla- 
tions and therefore the block maxima of the observables 
of the form g^-) asymptotically obeys one of the EVLs 
described above. We have already mentioned in the 
introduction that the behaviour of the system changes 
in the limit A — > A c as the correlations in the system 
increases [H[57]. Therefore, if the decay of correlations 
gets slower and slower the order of extreme statistics 



- in practise, the minimum length of our experimental 
bin from which we extract the maximum - needed to 
observe convergence to the predicted EVLs becomes 
higher and higher. In the limit, the bin length needed to 
decorrelate the data tends to infinity and this prevents 
from obtaining EVLs. At finite time, this effect gives rise 
to increasing deviations from the theoretical behaviour 
that can be explored numerically whenever the asymp- 
totic expected values for /i, a and k are known. This 
condition is easily met when dealing with systems whose 
physical measure is absolutely continuous with respect 
to Lebesgue. In this case the EVL parameters given in 
Eq. [6] only depend on the phase space dimension. For 
stochastically perturbed dynamical systems this is also 
true provided that we choose a suitable noise amplitude 
as explained in 119] . Moreover, once found the bin length 
needed to obtain the asymptotic results for A such that 
in the perturbed system we do not observe any critical 
transition, we can fix it to study the convergence to 
EVLs when approaching A c . 

In high dimensional systems, the theoretical frame- 
work described above is still valid even though it may 
be difficult to perform numerical simulations and control 
the variables in the phase space. In this case, critical 
transitions are highlighted by specific physical observ- 
ables that experience abrupt changes when the system 
crosses a tipping point. These observables undergo 
greater amplitude deviations in the direction of the 
state they are destined to shift to, as resulting from 
the bifurcation, than in the opposite direction, often 
showing an increase in the skewness of the distribution. 
Recently, the skewness increase near a tipping point has 
been proposed as an early warning indicator by Guttal 
et al. [28]. However, the method does not allow for 
having a quantitative estimation of the threshold value 
to be recognized as A c . An extreme value analysis on the 
maxima and on the minima of the distribution is capable 
of overcoming this problem and highlight the critical 
transition by providing a quantitative way to determine 
the value of A c . Let us call the physical observable </>, 
for A < A c , the extreme fluctuations of the observable 
are confined in the direction of the minima and of the 
maxima giving rise to a Weibull distribution typically 
observed for series bounded by an upper threshold. 
Instead, when A — > A c the fluctuations feel the presence 
of the other attractor in the direction of the state it is 
destined to shift to and this is reflected by a change 
in the extreme distributions from a Weibull one to a 
Gumbel one. This mechanism also allow to under- 
stand in which direction the transition will take place 
if the goal is to use the early warnings in a predictive way. 

The critical indicators introduced so far using the 
EVT results can be divided into two classes: dynami- 
cal indicators based on the divergence of the classical 
EVLs in the proximity of tipping points for the maxima 
of gi i = 1,2, 3 and physical indicators based on the 



5 



changes in the distributions of maxima (max(</>)) and "re- 
versed minima" (— min(<^>)) of a certain relevant physical 
observable <j> near the critical transitions. In the next sec- 
tion we will show an example in which the behaviour for 
the critical indicators introduced in this section can be 
derived analytically and observed numerically, defining 
guidelines for the application in other systems in Section 
5. 



IV. A CASE STUDY 

The system considered in the following analysis is 
widely used for modelling several natural phenomena fea- 
turing a bistable behaviour. 

The following stochastic differential equation (SDE): 

dx = -V'(x)dt + edW, (7) 

describes the over-damped motion of a material point in 
a potential under the effect of a stochastic forcing. The 
potential can be written as V(x) — \x A — ax 2 + cx with 
a, c > and it is represented in Fig. [IJfor some values of 
the parameters a and c. The stochastic forcing consists 
of a Wiener process W whose amplitude is modulated 
by the parameter e > 0. By letting e = 0, c = the 
system is completely deterministic and features two sta- 
ble fixed points (x~\ < and x~2 > ) and one unstable 
fixed point. It is trivial to check that extremes of any ob- 
servables extracted from the deterministic dynamics do 
not obey any of the classical EVLs. Instead, as soon as 
noise is switched on by setting a non-zero stochastic forc- 
ing, the physical measure becomes absolutely continuous 
with respect to Lebesgue so that, starting from any ini- 
tial conditions Xq, asymptotic EVL exists at any point £■ 
In fact, the probability distribution for the variable x can 
be computed by solving the Fokker Planck equation and 
considering a steady solution [29]. The result is obeyed 
in general and it is known as the Boltzman factor: 

^(^) = C( e )exp|-^M| (8) 

The dynamical behaviour of Eq. [7] can be understood 
by computing the mean exit time for the basin of attrac- 
tion of x~2 > as in [3D] , which for small values of c can 
be written as: 



(t(x 2 )) = 




while the mean exit time from the basin of attraction of 
Xi is obtained by changing c to — c. 
For c = the mean exit time from the negative and 
positive basin are the same. Instead for c > , fixing the 
initial condition in xq = £2, we notice that the mean exit 



time decreases when ,„ , c , /7 — s- 1. Once e is fixed such 
that < e <C a, the average exit time can be controlled 
by modulating c. 

Let us now analyse the behaviour of the dynamical 
indicators described in the previous section when c 
is changed. When noise amplitude is small and a 
local equilibrium approximation is used, the indicators 
converge towards the EVL as described in Eq. [6] with 
8(Q = 1 for all £ since it has been proven for the 
double well symmetric potential that the correlation 
decay is a stretched exponential (see Eqs. 47,48 in 
[31] ). As pointed out in the previous section, such a 
decay of correlations is enough to prove that gi extremes 
have classical EVLs. Instead, when (t(x 2 )) — > 0, the 
correlation decay becomes slower and slower (compare 
Fig. 4 and Fig. 5 in [31]) and the extremes will 
still obey a classical EVL but it becomes practically 
impossible to observe in a finite time series: this points 
at the importance of looking not only at the asymptotic 
statistics of processes, but also at the rate of convergence. 

As physical indicator we simply examine the maxima 
and reversed minima distribution of <p(x) — x which 
here represents the position of the particle in the right 
well of the potential V(X) being the initial condition 
xq = x~2- We remind that in a general case the analysis 
can be carried out by analysing the series of any relevant 
observables of the systems whose expectation value 
is different in the two basins of attraction. When 
c — > 0, e <C a, an harmonic expansion of the potential 
around x~2 holds and the particle is confined in both 
the direction of the minima and the maxima. In terms 
of extremes distributions we observe in both cases a 
Weibull type for the reasons described in the previous 
section. For increasing c the harmonic approximation 
of the potential does not hold any more and we have 
to consider the full expression of V(x): the maxima 
are still bounded by the quartic term but, towards the 
minima, the particle "feels" the presence of the other 
well experiencing a change of the distribution from 
Weibull to a Gumbel type. 



A. Numerical experiments 

The numerical experiments, performed on the one 
dimensional SDE described in Eq. [7j are devised to 
follow exactly the theoretical set-up described in this 
section. An ensemble of 100 realisations of the system 
is produced starting from x$ = x~2(c = 0), the noise 
amplitude is set at e = 0.5, dt = 0.1 under the potential 
V(x) — l/4x 4 — l/2x 2 + cx represented in Fig. [l] For 
each realisation, n = 1000 extremes are selected with 
two different bin length: m — 1000, m = 2000 for a total 
number of iterations s = 10 6 and s = 2 • 10 6 respectively. 
Even though for m = 1000, at c = 0, we get the shape 
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parameters as predicted from the theory, we repeat 
the experiment for increased s in order to show that 
the divergence at c cr a is not removed by considering 
longer time series. Approaching the critical transition, 
whenever the particle falls in the left well we interrupt 
the realisation reinitialising the initial condition xq 
counting, among all the realisations, the number of 
transitions experienced, which provides an indication 
on how likely a transition will occur during the time 
interval considered. Even if we expect our indicators 
to be effective as soon as transitions are observed, the 
reinitialisation helps in studying the behaviour of the in- 
dicators in the interval of c values for which the jump to 
the other basin of attraction is probable but not certain. 
This may help whenever only a realisation of the system 
is considered or when dealing with experimental data 
set. The inference procedure follows J32[ [33] where the 
authors have used the L-moments estimation described 
in [SI]. 

Let us first analyse what happens for the dynamical 
indicators for which results are shown in Fig. [2] where 
the shape parameters k for the gi observable are dis- 
played against the values of c. The left plots correspond 
to the case m = 10 3 , (3 = 1/3, the right ones to the case 
to = 10 4 , (3 = 1/3. Substituting the value of the param- 
eters chosen for the simulation in Eq[9j the exponential 
term vanishes when c — > 0.25 = c c ,it. Nonetheless, the 
value of c for which the system experiences at least one 
transition from one basin of attraction to another and 
a reinitialisation from the initial condition is needed is 
about c cr it — 0.21 in both the experiments. The differ- 
ence between the theoretical and experimental value of 
c cr it is due to the fact that the noise has a finite size. 
The divergence is far more evident for the observable g\ 
and 53 whereas it seems slower for the g 2 - This is due to 
the different weight assigned by g 2 to the minimum dis- 
tances: in fact, with respect to the observables g\ and <?3, 
g 2 weights more the points which come closer to £. Since 
spurious extremes are located relatively far away from 
they do not contribute sensibly to the divergence of 
the shape parameter of g 2 - This asymmetry with respect 
to the other observables may be removed by choosing a 
smaller /3 so that farther extremes are weighted more. 
Repeating the experiment setting j3 = 1/10 for g 2 we 
obtain an estimation of c cr u consistent with the other 
observables. To support this claim we have carried out 
the same experiments - not shown here - for the potential 
V(x) = ax 2 with a > 0, computing the statistics with re- 
spect to the point ( = 0. In this case, by setting a — >• 0, 
the potential gets flatter so that the material point can 
explore a wider region of the phase space causing the ap- 
pearance of spurious maxima in the distribution for finite 
size of m. This is reflected in the behaviour of g 2 which 
follows exactly the case of the double well potential so 
that for j3 — 1/3 the shape parameters diverges slower 
from the expected value than in the case f3 = 1/10. 
The divergent behaviour appears at lower values of c in 



the case to = 1000 with respect to the other case. In 
fact, by increasing the bin length the problem of lower 
decay of correlations appears at higher c than in the 
case to = 1000 as we improve the convergence by se- 
lecting more "authentic" extremes. Encouragingly, also 
in the to = 2000 experiment, the critical transition can 
be highlighted well before the probability of jumping to 
the other attracting states is higher than the 20%. From 
this example it is clear that is impossible to univocally 
define a threshold that indicates the tipping point. This 
is an intrinsic problem of the critical transitions in sys- 
tems driven by noise and it is directly related to the time 
scales on which we observe our system. Therefore the di- 
vergence obtained by the use of different gi should be crit- 
ically analysed to detect the approaching tipping point. 
Instead, in a modelling framework, one can test the sys- 
tem's behaviour with changing levels of noise. As dis- 
cussed in |21j when looking at average transitions rates, 
in a non one dimensional system, one should expect the 
same scaling properties of c cr u and e is the same as in 
the simple one dimensional model analysed here. 
The results for the physical indicators are shown in Fig[3] 
The upper plots represent the shape parameter for the 
maxima of <fi(x) = x, the middle plot the same but for 
the reversed minima whereas the lower plot is the same 
of Fig. 1 repeated here for convenience. The set-up is 
exactly the same as for the dynamical indicators. In this 
case what we observe is a change in the type of distri- 
bution only in the direction of the minima. This is ex- 
plainable by observing that an increase the bin length, 
the probability of observing a very low minimum within a 
single bin increases and this reflects in the overall distri- 
bution. In the case to = 2000 we observe even a Frechet 
law but for the values of the system for which the sys- 
tem has experienced a critical shift at least in the 20% 
of the realisations. In this case it seems indeed possible 
to define c cr u as the value for which the GEV distribu- 
tion of the minima or of the maxima changes sign. For 
more CPU demanding experiments or for experimental 
datasets an ensemble of realisations is usually not avail- 
able. In these cases the switching between different types 
of distribution for only the minima - or the maxima - may 
be still interpreted as a signal of an approaching tipping 
point using confidence intervals as uncertainty for the pa- 
rameters instead of the standard deviation of the sample. 
This approach we also applied to the analysis of turbu- 
lent energy data in the plane Couette flow about which 
we will report elsewhere. 



V. FINAL REMARKS 

This paper has addressed the problem of using 
extremes statistics to define robust indicators for ap- 
proaching tipping points of dynamical systems. The 
indicators have been grouped into two categories: 
dynamical indicators and physical indicators. For the 
former, by knowing some properties of the physical 
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measure it has been possible to identify the asymptotic 
EVLs when the system is well away from bifurcation 
points and thereby detect approaching tipping points 
by the divergence from the expected EVLs. Similarly, 
for the physical indicators constructed using extremes 
of relevant observables of the system, a comparison 
between maxima and reversed minima asymptotic 
EVL parameters provide a straightforward way to 
identify critical transitions. From the numerical results, 
described in the previous section, we suggest some 
guidelines to implement the indicators presented in an 
algorithmic way as described in Fig|4j 

• The dynamical indicators are mainly devised for 
applications in low dimensional systems as the tra- 
jectories must be explicitly computed. For de- 
terministic dynamical systems the local dimension 
should be known at the point £ whereas for stochas- 
tic dynamical systems the asymptotic EVLs depend 
only on the phase space dimension. The algorithm, 
represented schematically in the upper panel of Fig. 
|1J begins by setting A = and the computation of 
the length of the series s needed to obtain theo- 
retical parameters consistent with the asymptotic 
EVLs. Once s is determined, A is increased and the 
fitting procedure repeated until the GEV parame- 
ters diverge from the theoretical expected one. The 
critical A may be recognized when the experimental 
parameters are not consistent with the theoretical 
ones. 

• The physical indicators algorithm is described in 
the bottom panel of figure [4] It can be applied to a 
series of observables <p(x) originating from dynami- 
cal systems or to an experimental dataset. The first 
step, as in the previous case, is to chose a suitable 
s in order to obtain a Weibull law for maxima and 
reversed minima of <fr(x) in the unperturbed sys- 
tem. Once the length of the data series s is fixed, 
the critical A is the one for which the shape param- 
eter of the distribution either of the minima or of 
the maxima changes sign - provided that the other 
remains negative. 

The extreme values indicators present some noticeable 
advantages with respect to the indicators commonly 
used to highlight critical transitions: the knowledge of 
the asymptotic EVLs provides itself a robust indication 
of the order statistics needed to profit from an extreme 
value analysis. In other terms, we have the exact statis- 
tical model to conform to (under given conditions on the 
dynamical system's properties), and we can correctly 
infer that if the statistics does not obey the GEV model, 
then the conditions are not obeyed. This gives a much 
stronger mathematical framework to our analysis than 
in most previous investigations. In many applications 
this information is not available and this causes the 
identification of early warning signal as false tipping 



points. For the physical indicators, the interesting 
intuition that the skewness of the distribution changes 
when approaching the transition |28j can be related to 
modifications in the EVLs quantifiable with a change of 
sign of the shape parameters corresponding to a different 
type of extreme value distribution. Clearly, results will 
improve by choosing observables whose values change 
more in percentage crossing the tipping point. The strict 
connection between the dynamical indicators and the re- 
currences in a point suggest that we are able to highlight 
critical transitions only by looking at the behaviour of 
the systems in the neighbour of a specific point £ which 
can be located in any point of the attractor. This could 
help in all the situations for which the dynamic may 
be better represented or analysed in a sub-domain of 
the phase space. Moreover these indicators themselves 
provide interesting information on events which happen 
with very small probability but that are usually relevant 
in ecosystems, climate and financial models. Certainly 
the methods based on extreme value statistics require 
a great availability of good quality data and/or the 
possibility to perform time consuming simulations 
capable to extract authentic extremes. This problem is 
especially relevant for applications in climate science and 
in finance since the models involved are very demanding 
in terms of CPU time and long-term data availability 
is poor. However, such data will be more and more 
available in the future and not only for simplified models. 

The numerical tests carried out demonstrate the ap- 
plicability of the proposed indicators in practical appli- 
cations. Although the results meet the theoretical set- 
up and provide warning signals for approaching tipping 
points, it is indeed evident that other tests should be car- 
ried out to assess the general validity of this approach. 
We have already successfully tested the indicators on bi- 
dimensional systems obtaining results comparable to the 
ones shown in this paper and that will be reported else- 
where. Since the statistical tools used to study extreme 
values arc commonly distributed with scientific software, 
the algorithm can be easily checked and compared to 
other methods. Our aim is to test the indicator on com- 
plex systems arising from fluid dynamical studies which 
feature tipping points whose nature remain not well un- 
derstood and on other theoretical low dimensional mod- 
els. One must bear in mind that, as discussed in [2T] . 
great caution must be paid when trying to extend con- 
siderations valid for one dimensional models onto higher 
dimensional systems: the operation of defining a one di- 
mensional effective projected dynamics is far from being 
trivial. 
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FIG. 1. Potential V(x) = l/4x 4 - 



l/2x 2 + cx for different values of parameter c. 
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FIG. 2. Extreme value shape parameter k vs the potential asymmetry parameter c for an ensemble of 30 trajectories starting 
in xo = X2 of the system described in Eq. [7] (a): g\ observable, (b): g2 observable with /3 = 1/3, (c): g^ observable with 
/3 = 1/3, (7 = 1, (d): percentage of runs for which the particle jump at least once in the left% well. Left: n = 1000, m = 1000. 
Right: n = 1000, m — 2000. The errorbars represent a standard deviation of the sample, the red bars represent theoretical 
expected parameters for c = case. 
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FIG. 3. Extreme value shape parameter k vs the potential asymmetry parameter c for an ensemble of 30 trajectories starting 
in :ro = X2 of the system described in Eq. [7| (a) : maxima for the observable 4>{x) = x (b) : reversed minima for the observable 
4>{x) — x\ (c): percentage of runs for which the particle jump at least once in the left% well. Left: n = 1000, m = 1000. 
Right: n = 1000, m = 2000. The errorbars represent a standard deviation of the sample, the red bars represent the Gumbel 
distribution (k = 0). 
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Dynamical system: Xj=f J (x ,A) 



1) Choose an initial condition x Q . 

2) Produce the series of x , x r ..,x-,... 

3) Choose a point Z, on the attractor. 




Block Maxima approach 



1) Compute the series X-=g(dlst(f i (x ,M,O) ■ 

2) Divide the series X Q ,X 1 X s into n bins 

each containing m observations. 

3) Take maxima of X- in each bin. 

4) Fit the maxima to a GEV model to get k. 





Dynamical system: x^f^x^A) 



1) Choose an initial condition x Q 

2) Produce the series of x Q , x 1 ...,x-,..., 

3) Choose a point T, on the attractor. 



Block Maxima approach 



1) Compute the series X.=q(dist(f>(x ,>>),Q) . 

2) Divide the series X Q ,X X X s into n bins 

each containing m observations. 

3) Take maxima of X. in each bin. 

4) Fit the maxima to a GEV model to get k. 




Dynamical system: x-- 



-■f } (x ,A) 



1) Choose an initial condition x Q . 

2) Compute q>(x Q ), q>(x 1 )..., <p(x s ). 



Block Maxima approach 



1) Divide the series <p(x Q ), <p(Xj )..., <p(xj into n bins 
each containing m observations. 

2) take maxima and minima of <p(Xj ) in each bin. 

3) Fit the maxima to a GEV model to get k mgx . 

4) Fit the reversed minima to a GEV, get k ■ . 





Dynamical system: x-=f J (x ,A) 



1) Choose an initial condition x Q . 

2) Compute q>(x Q ), q>(x 2 )...,q>(x s ). 



Block Maxima approach 

1) Divide the series <p(Xq), (p(x 1 )..., (pfx^ into n bins 

each containing m observations. 

3) take maxima and minima of <p(x, ) in each bin. 

4) Fit the maxima to a GEV model to get k mgx . 

5) Fit the reversed minima to a GEV, get k . 



-H end 



FIG. 4. Schematic representation of the numerical procedure to be used to infer the critical transition happening at A = A c . 
Upper panel: dynamical indicators. Lower panel: physical indicators. See the text for further explanations 



